The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file
Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd
|
| | 0%
|
|....... | 4%
|
|.............. | 9%
|
|...................... | 13%
|
|............................. | 17%
|
|.................................... | 22%
|
|........................................... | 26%
|
|................................................... | 30%
|
|.......................................................... | 35%
|
|................................................................. | 39%
|
|........................................................................ | 43%
|
|............................................................................... | 48%
|
|....................................................................................... | 52%
|
|.............................................................................................. | 57%
|
|..................................................................................................... | 61%
|
|............................................................................................................ | 65%
|
|................................................................................................................... | 70%
|
|........................................................................................................................... | 74%
|
|.................................................................................................................................. | 78%
|
|......................................................................................................................................... | 83%
|
|................................................................................................................................................ | 87%
|
|........................................................................................................................................................ | 91%
|
|............................................................................................................................................................... | 96%
|
|......................................................................................................................................................................| 100%
output file: /tmp/RtmphhTqE3/file42a957179148
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Attaching package: ‘ftExtra’
The following object is masked from ‘package:flextable’:
separate_header
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes_Remake %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
seqDep <- countMat %>% apply(1, sum)
min(seqDep)
[1] 852
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20
4.377867 5.122635 4.697075 5.919422 5.091001 3.868273 5.527489 4.512680 4.845492 4.782770 5.290467 4.835734 4.308318 4.664723 4.653949 5.022267
3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500
5.364248 4.956451 4.844607 3.744475 5.000444 4.749852 4.918364 5.142239 4.981104 4.204087 4.331178 4.789362 3.432382 5.736709 5.205365 5.265310
3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500
5.492566 4.997043 4.849876 4.869032 4.255819 4.360230 4.996300 4.423317 4.437111 4.771900 4.308965 4.327765 4.891794 4.683325 5.170850 4.084572
4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500
4.491632 2.805419 4.592134 4.741784 4.721501 4.533240 4.428191 4.456642 4.070497 3.924195 4.118489 4.097255 4.642190 5.118735 5.533426 5.073586
5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
4.951162 4.296940 4.977480 4.909048 4.280855 4.186882 4.900147 4.860817 5.399822 4.662039 5.209531
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2
0.011467055 0.007103701 0.012255238 0.003349956 0.008110994 0.041006425 0.007615813 0.011728862 0.009093947 0.008972655 0.008846283 0.007282732 0.012661180 0.009597579
3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2
0.011329508 0.005201462 0.008190891 0.009811517 0.010195269 0.021596351 0.008966390 0.010189900 0.008947565 0.008349931 0.009991783 0.014580123 0.012148164 0.008882062
3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5
0.063562627 0.005282231 0.009346957 0.006143886 0.005836946 0.010690993 0.008613137 0.009301306 0.010889362 0.012524899 0.007823094 0.011012647 0.008667124 0.009030647
4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180
0.005918521 0.010062620 0.009452742 0.014296051 0.006724997 0.014380807 0.011650840 0.124685269 0.011358859 0.008391334 0.010995971 0.011373619 0.011834659 0.012159030
5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2
0.014452214 0.011299615 0.014013080 0.013216952 0.009980342 0.008042424 0.005883545 0.010729090 0.008357545 0.011428033 0.011626609 0.008644633 0.015737093 0.008575283
C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
0.007366125 0.006551531 0.004008265 0.007772224 0.005160141
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied observed species numbers
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-219.518 -35.649 -1.039 49.226 201.182
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 248133.486 132566.735 1.872 0.0655 .
log(Size_Class) 33.680 8.032 4.193 8.02e-05 ***
I(log(Size_Class)^2) -6.473 1.494 -4.332 4.92e-05 ***
lat -12915.027 6907.881 -1.870 0.0658 .
I(lat^2) 168.133 89.927 1.870 0.0658 .
depth 3.905 3.181 1.227 0.2238
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 75.51 on 69 degrees of freedom
Multiple R-squared: 0.2656, Adjusted R-squared: 0.2124
F-statistic: 4.991 on 5 and 69 DF, p-value: 0.0005849
Rarified chao1 estimates
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-521.02 -131.00 -37.98 109.00 1052.27
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 795248.56 436480.86 1.822 0.072796 .
log(Size_Class) 89.88 26.45 3.399 0.001128 **
I(log(Size_Class)^2) -17.77 4.92 -3.611 0.000574 ***
lat -41431.39 22744.45 -1.822 0.072850 .
I(lat^2) 539.62 296.09 1.822 0.072715 .
depth 16.95 10.47 1.618 0.110139
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 248.6 on 69 degrees of freedom
Multiple R-squared: 0.1986, Adjusted R-squared: 0.1405
F-statistic: 3.419 on 5 and 69 DF, p-value: 0.008097
As above but without latitude and depth
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-468.58 -139.15 -22.35 112.75 1118.50
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 539.262 44.634 12.082 < 2e-16 ***
log(Size_Class) 90.573 26.507 3.417 0.001045 **
I(log(Size_Class)^2) -18.060 4.931 -3.663 0.000473 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 249.5 on 72 degrees of freedom
Multiple R-squared: 0.1577, Adjusted R-squared: 0.1343
F-statistic: 6.74 on 2 and 72 DF, p-value: 0.002074
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.32735 -0.16629 0.02731 0.30484 0.76295
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.419e+03 7.933e+02 1.788 0.0782 .
log(Size_Class) 2.055e-01 4.807e-02 4.276 6.00e-05 ***
I(log(Size_Class)^2) -3.699e-02 8.943e-03 -4.136 9.79e-05 ***
lat -7.362e+01 4.134e+01 -1.781 0.0793 .
I(lat^2) 9.580e-01 5.382e-01 1.780 0.0795 .
depth 1.321e-02 1.904e-02 0.694 0.4899
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4519 on 69 degrees of freedom
Multiple R-squared: 0.2905, Adjusted R-squared: 0.239
F-statistic: 5.649 on 5 and 69 DF, p-value: 0.0002011
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.014081 -0.004516 -0.001758 0.000707 0.099175
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.476e+01 2.619e+01 -0.946 0.3476
log(Size_Class) -4.146e-03 1.587e-03 -2.613 0.0110 *
I(log(Size_Class)^2) 7.087e-04 2.952e-04 2.401 0.0191 *
lat 1.289e+00 1.365e+00 0.945 0.3480
I(lat^2) -1.676e-02 1.776e-02 -0.944 0.3486
depth -3.746e-04 6.284e-04 -0.596 0.5531
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01492 on 69 degrees of freedom
Multiple R-squared: 0.1063, Adjusted R-squared: 0.04159
F-statistic: 1.642 on 5 and 69 DF, p-value: 0.1605
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
219.0394 219.0394 191.6764 191.6764 200.2007 157.3011 157.3011 183.0517 202.9587 295.9379 295.9379 268.5750 268.5750 277.0992 234.1997 234.1997 259.9502 259.9502 327.4518
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
327.4518 300.0888 300.0888 308.6131 265.7136 265.7136 291.4641 311.3712 311.3712 332.8183 332.8183 305.4553 305.4553 313.9796 313.9796 271.0800 271.0800 296.8306 296.8306
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
321.6979 294.3349 294.3349 302.8592 302.8592 259.9597 259.9597 259.9597 285.7102 285.7102 305.6173 305.6173 290.3570 290.3570 262.9941 262.9941 271.5184 271.5184 228.6188
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
228.6188 228.6188 254.3694 254.3694 274.2764 274.2764 249.3269 221.9639 221.9639 230.4882 230.4882 187.5886 187.5886 187.5886 213.3392 213.3392 233.2462 233.2462
$se.fit
[1] 26.41573 26.41573 25.69552 25.69552 25.55840 27.19115 27.19115 28.71794 33.31516 18.53028 18.53028 17.92699 17.92699 16.83554 19.61469 19.61469 21.04649 21.04649 18.84519
[20] 18.84519 18.39317 18.39317 16.87013 19.66040 19.66040 20.75488 27.57872 27.57872 18.84306 18.84306 18.36078 18.36078 16.60167 16.60167 19.19454 19.19454 20.13853 20.13853
[39] 18.16789 17.54242 17.54242 15.65854 15.65854 18.05154 18.05154 18.05154 19.02151 19.02151 25.99078 25.99078 18.90217 18.90217 18.03402 18.03402 16.32633 16.32633 18.03550
[58] 18.03550 18.03550 19.08494 19.08494 25.56051 25.56051 23.94049 23.00762 23.00762 21.87070 21.87070 22.64968 22.64968 22.64968 23.62789 23.62789 28.60636 28.60636
$df
[1] 69
$residual.scale
[1] 75.50836
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
414.2733 414.2733 319.9690 319.9690 396.7495 263.9882 263.9882 373.9864 352.6608 620.7464 620.7464 526.4421 526.4421 603.2225 470.4612 470.4612 580.4595 580.4595 703.5800
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
703.5800 609.2758 609.2758 686.0562 553.2949 553.2949 663.2932 641.9675 641.9675 714.7463 714.7463 620.4421 620.4421 697.2225 697.2225 564.4612 564.4612 674.4595 674.4595
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
681.7164 587.4121 587.4121 664.1925 664.1925 531.4312 531.4312 531.4312 641.4295 641.4295 620.1038 620.1038 592.5451 592.5451 498.2408 498.2408 575.0212 575.0212 442.2600
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
442.2600 442.2600 552.2582 552.2582 530.9325 530.9325 477.2943 382.9901 382.9901 459.7705 459.7705 327.0092 327.0092 327.0092 437.0074 437.0074 415.6818 415.6818
$se.fit
[1] 86.97476 86.97476 84.60346 84.60346 84.15196 89.52786 89.52786 94.55488 109.69140 61.01161 61.01161 59.02528 59.02528 55.43162 64.58209 64.58209 69.29635
[18] 69.29635 62.04847 62.04847 60.56018 60.56018 55.54551 64.73259 64.73259 68.33621 90.80394 90.80394 62.04146 62.04146 60.45353 60.45353 54.66162 54.66162
[35] 63.19872 63.19872 66.30685 66.30685 59.81845 57.75905 57.75905 51.55631 51.55631 59.43537 59.43537 59.43537 62.62903 62.62903 85.57562 85.57562 62.23608
[52] 62.23608 59.37768 59.37768 53.75503 53.75503 59.38253 59.38253 59.38253 62.83786 62.83786 84.15892 84.15892 78.82495 75.75342 75.75342 72.01009 72.01009
[69] 74.57490 74.57490 74.57490 77.79571 77.79571 94.18749 94.18749
$df
[1] 69
$residual.scale
[1] 248.6141
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
4.460000 4.460000 4.316870 4.316870 4.263512 4.006600 4.006600 4.096129 4.380305 4.922816 4.922816 4.779687 4.779687 4.726328 4.469417 4.469417 4.558945 4.558945 5.121512
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
5.121512 4.978383 4.978383 4.925024 4.668113 4.668113 4.757641 5.041818 5.041818 5.170247 5.170247 5.027117 5.027117 4.973759 4.973759 4.716847 4.716847 4.806376 4.806376
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
5.119395 4.976266 4.976266 4.922907 4.922907 4.665996 4.665996 4.665996 4.755524 4.755524 5.039701 5.039701 4.956219 4.956219 4.813089 4.813089 4.759731 4.759731 4.502820
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
4.502820 4.502820 4.592348 4.592348 4.876525 4.876525 4.735050 4.591920 4.591920 4.538562 4.538562 4.281650 4.281650 4.281650 4.371179 4.371179 4.655356 4.655356
$se.fit
[1] 0.15808462 0.15808462 0.15377457 0.15377457 0.15295393 0.16272512 0.16272512 0.17186218 0.19937421 0.11089422 0.11089422 0.10728388 0.10728388 0.10075207 0.11738388
[16] 0.11738388 0.12595249 0.12595249 0.11277880 0.11277880 0.11007370 0.11007370 0.10095908 0.11765743 0.11765743 0.12420735 0.16504451 0.16504451 0.11276605 0.11276605
[31] 0.10987985 0.10987985 0.09935253 0.09935253 0.11486948 0.11486948 0.12051879 0.12051879 0.10872554 0.10498238 0.10498238 0.09370833 0.09370833 0.10802925 0.10802925
[46] 0.10802925 0.11383402 0.11383402 0.15554155 0.15554155 0.11311979 0.11311979 0.10792438 0.10792438 0.09770471 0.09770471 0.10793321 0.10793321 0.10793321 0.11421358
[61] 0.11421358 0.15296658 0.15296658 0.14327160 0.13768881 0.13768881 0.13088496 0.13088496 0.13554674 0.13554674 0.13554674 0.14140086 0.14140086 0.17119443 0.17119443
$df
[1] 69
$residual.scale
[1] 0.4518789
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.020303322 0.020303322 0.022872940 0.022872940 0.021919895 0.025509948 0.025509948 0.022656554 0.019388204 0.011062673 0.011062673 0.013632292 0.013632292 0.012679247
15 16 17 18 19 20 21 22 23 24 25 26 27 28
0.016269300 0.016269300 0.013415905 0.013415905 0.006958323 0.006958323 0.009527941 0.009527941 0.008574897 0.012164950 0.012164950 0.009311555 0.006043205 0.006043205
29 30 31 32 33 34 35 36 37 38 39 40 41 42
0.005735595 0.005735595 0.008305213 0.008305213 0.007352168 0.007352168 0.010942221 0.010942221 0.008088827 0.008088827 0.006506674 0.009076292 0.009076292 0.008123248
43 44 45 46 47 48 49 50 51 52 53 54 55 56
0.008123248 0.011713300 0.011713300 0.011713300 0.008859906 0.008859906 0.005591556 0.005591556 0.009378051 0.009378051 0.011947669 0.011947669 0.010994624 0.010994624
57 58 59 60 61 62 63 64 65 66 67 68 69 70
0.014584677 0.014584677 0.014584677 0.011731283 0.011731283 0.008462932 0.008462932 0.013402419 0.015972038 0.015972038 0.015018993 0.015018993 0.018609046 0.018609046
71 72 73 74 75
0.018609046 0.015755652 0.015755652 0.012487301 0.012487301
$se.fit
[1] 0.005218045 0.005218045 0.005075779 0.005075779 0.005048692 0.005371218 0.005371218 0.005672814 0.006580929 0.003660388 0.003660388 0.003541218 0.003541218 0.003325617
[15] 0.003874598 0.003874598 0.004157430 0.004157430 0.003722594 0.003722594 0.003633304 0.003633304 0.003332450 0.003883628 0.003883628 0.004099827 0.005447777 0.005447777
[29] 0.003722173 0.003722173 0.003626906 0.003626906 0.003279421 0.003279421 0.003791603 0.003791603 0.003978075 0.003978075 0.003588804 0.003465251 0.003465251 0.003093117
[43] 0.003093117 0.003565821 0.003565821 0.003565821 0.003757425 0.003757425 0.005134104 0.005134104 0.003733850 0.003733850 0.003562360 0.003562360 0.003225030 0.003225030
[57] 0.003562651 0.003562651 0.003562651 0.003769953 0.003769953 0.005049109 0.005049109 0.004729098 0.004544822 0.004544822 0.004320241 0.004320241 0.004474116 0.004474116
[71] 0.004474116 0.004667349 0.004667349 0.005650773 0.005650773
$df
[1] 69
$residual.scale
[1] 0.01491559
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
bold(i = ~ p< 0.05, j = "p") %>%
colformat_md() %>%
set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.5 x 105 | 1.3 x 105 | 1.87 | 0.065 |
log(Size Class) | 3.4 x 101 | 8.0 x 100 | 4.19 | < 0.001 | |
log(Size Class)2 | -6.5 x 100 | 1.5 x 100 | -4.33 | < 0.001 | |
Latitude | -1.3 x 104 | 6.9 x 103 | -1.87 | 0.066 | |
Latitude2 | 1.7 x 102 | 9.0 x 101 | 1.87 | 0.066 | |
Depth | 3.9 x 100 | 3.2 x 100 | 1.23 | 0.224 | |
Richness (Chao1) | Intercept | 8.0 x 105 | 4.4 x 105 | 1.82 | 0.073 |
log(Size Class) | 9.0 x 101 | 2.6 x 101 | 3.40 | 0.001 | |
log(Size Class)2 | -1.8 x 101 | 4.9 x 100 | -3.61 | < 0.001 | |
Latitude | -4.1 x 104 | 2.3 x 104 | -1.82 | 0.073 | |
Latitude2 | 5.4 x 102 | 3.0 x 102 | 1.82 | 0.073 | |
Depth | 1.7 x 101 | 1.0 x 101 | 1.62 | 0.110 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.9 x 102 | 1.79 | 0.078 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.28 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.14 | < 0.001 | |
Latitude | -7.4 x 101 | 4.1 x 101 | -1.78 | 0.079 | |
Latitude2 | 9.6 x 10-1 | 5.4 x 10-1 | 1.78 | 0.079 | |
Depth | 1.3 x 10-2 | 1.9 x 10-2 | 0.69 | 0.490 | |
Evenness (Pielou J) | Intercept | -2.5 x 101 | 2.6 x 101 | -0.95 | 0.348 |
log(Size Class) | -4.1 x 10-3 | 1.6 x 10-3 | -2.61 | 0.011 | |
log(Size Class)2 | 7.1 x 10-4 | 3.0 x 10-4 | 2.40 | 0.019 | |
Latitude | 1.3 x 100 | 1.4 x 100 | 0.94 | 0.348 | |
Latitude2 | -1.7 x 10-2 | 1.8 x 10-2 | -0.94 | 0.349 | |
Depth | -3.7 x 10-4 | 6.3 x 10-4 | -0.60 | 0.553 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013935 -0.005053 -0.002494 0.000907 0.105945
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.980e+01 2.751e+01 -0.720 0.4743
log(Size_Class) -3.292e-03 1.667e-03 -1.975 0.0523 .
I(log(Size_Class)^2) 5.754e-04 3.102e-04 1.855 0.0679 .
lat 1.030e+00 1.434e+00 0.718 0.4751
I(lat^2) -1.338e-02 1.866e-02 -0.717 0.4760
depth -2.432e-04 6.603e-04 -0.368 0.7138
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01567 on 69 degrees of freedom
Multiple R-squared: 0.06742, Adjusted R-squared: -0.0001586
F-statistic: 0.9977 on 5 and 69 DF, p-value: 0.4258
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
0.0116609352 0.0116609352 0.0135814780 0.0135814780 0.0133501419 0.0160295241 0.0160295241 0.0139649206 0.0099302928 0.0042917835 0.0042917835 0.0062123262
13 14 15 16 17 18 19 20 21 22 23 24
0.0062123262 0.0059809902 0.0086603723 0.0086603723 0.0065957688 0.0065957688 0.0010653312 0.0010653312 0.0029858739 0.0029858739 0.0027545378 0.0054339200
25 26 27 28 29 30 31 32 33 34 35 36
0.0054339200 0.0033693165 -0.0006653113 -0.0006653113 0.0001751617 0.0001751617 0.0020957044 0.0020957044 0.0018643683 0.0018643683 0.0045437505 0.0045437505
37 38 39 40 41 42 43 44 45 46 47 48
0.0024791470 0.0024791470 0.0008731389 0.0027936816 0.0027936816 0.0025623455 0.0025623455 0.0052417277 0.0052417277 0.0052417277 0.0031771242 0.0031771242
49 50 51 52 53 54 55 56 57 58 59 60
-0.0008575036 -0.0008575036 0.0032944705 0.0032944705 0.0052150133 0.0052150133 0.0049836772 0.0049836772 0.0076630594 0.0076630594 0.0076630594 0.0055984559
61 62 63 64 65 66 67 68 69 70 71 72
0.0055984559 0.0015638281 0.0015638281 0.0066369468 0.0085574895 0.0085574895 0.0083261534 0.0083261534 0.0110055356 0.0110055356 0.0110055356 0.0089409321
73 74 75
0.0089409321 0.0049063043 0.0049063043
$se.fit
[1] 0.005482717 0.005482717 0.005333235 0.005333235 0.005304774 0.005643660 0.005643660 0.005960553 0.006914730 0.003846052 0.003846052 0.003720837 0.003720837 0.003494300
[15] 0.004071127 0.004071127 0.004368305 0.004368305 0.003911413 0.003911413 0.003817594 0.003817594 0.003501480 0.004080615 0.004080615 0.004307780 0.005724101 0.005724101
[29] 0.003910971 0.003910971 0.003810871 0.003810871 0.003445761 0.003445761 0.003983922 0.003983922 0.004179853 0.004179853 0.003770837 0.003641017 0.003641017 0.003250008
[43] 0.003250008 0.003746688 0.003746688 0.003746688 0.003948011 0.003948011 0.005394518 0.005394518 0.003923240 0.003923240 0.003743051 0.003743051 0.003388611 0.003388611
[57] 0.003743358 0.003743358 0.003743358 0.003961175 0.003961175 0.005305212 0.005305212 0.004968969 0.004775346 0.004775346 0.004539374 0.004539374 0.004701055 0.004701055
[71] 0.004701055 0.004904088 0.004904088 0.005937394 0.005937394
$df
[1] 69
$residual.scale
[1] 0.01567214
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.9 x 102 | 1.79 | 0.078 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.28 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 8.9 x 10-3 | -4.14 | < 0.001 | |
Latitude | -7.4 x 101 | 4.1 x 101 | -1.78 | 0.079 | |
Latitude2 | 9.6 x 10-1 | 5.4 x 10-1 | 1.78 | 0.079 | |
Depth | 1.3 x 10-2 | 1.9 x 10-2 | 0.69 | 0.490 | |
Evenness (Pielou J) | Intercept | -2.0 x 101 | 2.8 x 101 | -0.72 | 0.474 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.052 | |
log(Size Class)2 | 5.8 x 10-4 | 3.1 x 10-4 | 1.86 | 0.068 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.72 | 0.475 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.72 | 0.476 | |
Depth | -2.4 x 10-4 | 6.6 x 10-4 | -0.37 | 0.714 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)